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We investigate the structure and dynamics of sodium disihcate by 
means of molecular dynamics computer simulation. We show that 
the structure is described by a partially destroyed tetrahedral Si04 
network and a spherical super structure formed by the silicon and 
sodium atoms. The static structure factor of our simulation is in very 
good agreement with one from a neutron scattering experiment. For 
1008 particles we find strong finite size effects in the dynamics which 
are due to the missing of modes contributing to the boson peak. 



1 Introduction 

In recent years several molecular dynamics computer simulations have been done 
in order to investigate the structure and dynamics of sodium silicate melts and 
glasses (Smith, Greaves, and Gillan 1995, Cormack and Cao 1997). By using the 
potential proposed by Vessal, Amini, Fincham and Catlow (1989) these authors 
have found that the structure of systems like, e.g., sodium disilicate (SDS) is 
characterized by a microsegregation in which the sodium atoms form clusters of 
a few atoms between bridged Si04 units. In order to see whether this somewhat 
surprising result is reproduced also by a different model from the one of Vessal 
et al. (1989) we have performed simulations of SDS using a different potential 
(discussed in detail below). In addition to the investigation of the structure we 
also study the dynamical properties of SDS in order to see whether the finite size 
effects that have been observed in pure silica (Horbach, Kob, Binder, and Angell 
1996, Horbach, Kob, and Binder 1999a, and Horbach, Kob, and Binder 1999b) are 
present in SDS as well. 

2 Model 

The model potential we use to describe the interactions between the ions in SDS 
is the one proposed by Kramer, de Man, and van Santen (1991) which is a gener- 
alization of the so called BKS potential (van Beest, Kramer, and van Santen 1990) 
for pure silica. It has the following functional form: 



+ A^pe-^'^^^-^ «,/5g [Si,Na,0]. (1) 



Here r is the distance between an ion of type a and an ion of type (3. The values 
of the parameters Aajj, Bap and Cap can be found in the original publication. The 
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potential (|1]) has been optimized by Kramer et al. for zeolites, i.e. for systems 
that have Al ions in addition to Si, Na and O. In that paper the authors used 
for silicon and oxygen the partial charges gsi = 2.4 and go = —1-2, respectively, 
whereas sodium was assigned its real ion charge gNa = 1-0. With this choice 
charge neutrality is not fulfilled in systems like SDS. To overcome this problem we 
introduced for the sodium ions a position dependent charge q{r) instead of gNa, 



q[r) 



0.6(l + ln[C(rc-r)Vl]) r < 
0.6 r > Tr 



which means that for r > Vc charge neutrality is valid (g(r) = 0.6 for r > Vc). Note 
that g(r) is continuous at Vc- We have fixed the parameters Vc and C such that the 
experimental mass density of SDS and the static structure factor from a neutron 
scattering experiment (see below) are reproduced well. From this fitting we have 
obtained the values = 4.9A and C = 0.0926A~^. With this choice the charge 
g(r) crosses smoothly over from g(r) = 1.0 at 1.7 A to g(r) = 0.6 for r > Tc. 

The simulations have been done at constant volume with the density of the 
system fixed to 2.37 g/cm^. The equations of motion were integrated with the 
velocity form of the Verlet algorithm and the Coulombic contributions to the po- 
tential and the forces were calculated via Ewald summation. The time step of the 
integration was 1.6 fs. In this paper we investigate the properties of SDS in the 
liquid state at T = 2100 K and in the glass state at T = 300 K. The equilibration 
time at T = 2100 K was two million time steps thus corresponding to a real time 
of 3.5 ns. At this temperature we simulated systems with N = 1008 and = 8064 
particles. In order to improve the statistics two independent runs were done for the 
large system and eight independent runs for the small system. The glass state was 
produced by cooling the system from equilibrated configurations at T = 1900 K 
with a cooling rate of 1.16 ■ 10^^ K/s. The pressure is 4.5 GPa at T = 2100 K and 
0.96 GPa at T = 300 K. 



3 Results 

In order to demonstrate that our model is able to reproduce the structure of real 
SDS very well we compare the static structure factor S'^'^^{q) with the one from a 
neutron scattering experiment by Misawa, Price, and Suzuki (1980). To calculate 
S^^^{q) one has to weight the partial structure factors from the simulation with 
the experimental coherent neutron scattering lengths ba (« G [Si, Na, O]): 

^°^"(?) = ^^-^ E (e^'^ I^-^'l) . (3) 

The values for are 0.4149 • lO^^^ pm, 0.363 ■ lO'^^ and 0.5803 ■ 10"^^ 
for silicon, sodium and oxygen, respectively. They are taken from Susman, Volin, 
Montague, and Price (1991) for silicon and oxygen and from Bacon (1972) for 
sodium. Fig. |I] shows 5'°'^"(g) from the simulation and the experiment at T = 
300 K. We see that the overall agreement between simulation and experiment is 
good. For g > 2.3 A~^, which corresponds to length scales of next nearest Si-0 
and Na-0 neighbors, the largest discrepancy is at the peak located at g = 2.8 



where the simulation underestimates the experiment by approximately 15% in 
amplitude. Very well reproduced is the peak at g = 1.7 A~^, which is called the 
first sharp diffraction peak, and which is a prominent feature in pure silica as 
well. In silica this peak arises from the tetrahedral network structure since the 
length scale which corresponds to it, i.e. 2tx/1.7 A = 3.7 A, is approximately 
the spatial extent of two connected Si04 tetrahedra. From the figure we recognize 
that this structure is partly present in SDS also. The peak at g = 0.95 is 
not present in the experimental data which might be due to the fact that in this q 
region the experimental resolution is not sufficient. By looking at the coordination 
number distributions, discussed below, we see that the peak at g = 0.95 is 
related to a super structure which is formed by the sodium and silicon atoms. In 
agreement with this interpretation the length scale corresponding to this peak, 
i.e. 27r/0.95 A = 6.6 A, is two times the mean distance of nearest Na-Na or 
Na-Si neighbors. 

The coordination number distribution Paplz) for different pairs ajS gives the 
probability that an ion of type a has exactly z nearest neighbors of type f3. By 
definition two neighboring atoms have a distance from each other which is less than 
the location of the first minimum rmin of the corresponding partial pair correlation 
function Qapir). From the functions Qapij) we find for Tmin the values 3.6 A, 
5.0 A, 2.35 A, 5.0 A, 3.1 A, and 3.15 A for the Si-Si, Si-Na, Si-0, Na-Na, Na- 
O, and 0-0 correlations. We note that Psi^o{z) is larger than 0.99 for z = A 
at T = 2100 K and T = 300 K which means that nearly every silicon atom is 
four fold coordinated with oxygen atoms forming a Si04 tetrahedron. Some of 
the distribution functions are shown in Fig. ^ We recognize from Fig. ^ that 
about 65% of the oxygen atoms are bridging oxygens between two tetrahedra 
(Po-Si(-2 = 2) ^ 0.65), and that about 28% of the oxygen atoms form dangling 
bonds {z = 1) with corresponding silicon atoms. In the neighborhood of these 
dangling bonds sodium atoms are located. This means that the sodium atoms 
partly destroy the (disordered) tetrahedral network which is the structure for pure 
silica. A significant number of oxygen atoms have no silicon atoms, but only 
sodium atoms as direct neighbors. In Fig. ^jb we show -PNa-Na(-2) at T = 2100 K 
and T = 300 K, and we see that essentially the two distributions coincide with a 
a mean value between z = 8 and z = 9. Basically the same distribution is found 
for P^a-si{z)- Therefore, every sodium atom is surrounded by 8-9 other sodium 
atoms and 8-9 silicon atoms. Since the mean distance between Na-Na and Na-Si 
neighbors is approximately the same, i.e. about 3.3 A, we can conclude that sodium 
and silicon atoms form a spherical super structure in which every sodium atom is 
surrounded by a first shell of oxygen atoms and a second shell of silicon and sodium 
atoms. In Fig. ||b we have also included PNa-o(-2^) and we recognize that every 
sodium atom has on average about 4-5 nearest oxygen neighbors. There is again 
no essential difference between Pne-o for T = 2100 K and Pne-o for T = 300 K, 
although the pressure is about a factor 4.5 higher at T = 2100 K. This means that 
the structural features which we observe at T = 2100 K are not formed due to 
the relatively high pressure. Nevertheless, we emphasize that the small difference 
between P{z) in the liquid and in the glass state is partly due to the high cooling 
rate of about 10^^ K/s we used to produce the structures at T = 300 K. A careful 
analysis of the cooling rate effects for SDS will be presented elsewhere (Horbach, 
Kob, and Binder 1999c). 



Having described the structure of SDS we turn now our attention to a dynam- 
ical quantity, namely the self part of the dynamic structure factor Ss{q, v)-, which 
depends on frequency z/ and the magnitude of the wave- vector q. It is defined by 

AT /"OO , . 

S^[q^y) = ^^g-2..t/g^q.[r.(t)-r.(0)]\ 

iV J —oo 
AT poo 

= -f / dte-'^'''F,{q,t) ae[Si,Na,0] (4) 

iV J — oo 

where Ta(t) is the position vector of a tagged particle of type a at time t and Fs{q, t) 
is the incoherent intermediate scattering function. The details of the Fourier trans- 
formation in (ID are given elsewhere (Horbach et al. 1999a). Ss{qiv) for silicon, 
sodium and oxygen is shown in Fig. |^a at the temperature T = 2100 K and the 
wave-vector q = 1.7 for the two system sizes N = 1008 and = 8064. 
For sodium the curves for the small and the large system coincide over the whole 
frequency range. This is not the case for silicon and oxygen for which Ssiq,!^) 
has no system size dependence for frequencies u > 1.1 THz whereas for smaller 
frequencies there is a missing of intensity for the small system. The vibrational 
modes causing a shoulder around z/ = 0.9 THz in Ss{q, i') are usually called boson 
peak excitations. In the small system a part of these excitations is missing due 
to the loss of intensity for v < 1.1 THz. An explanation of this behavior for the 
case of silica, which shows qualitatively the same finite size effects, can be found 
in Horbach et al. (1999a and 1999b). Due to the sum rule / dv Ss{q,i^) = 1 the 
missing of the boson peak excitations for frequencies 0.4 THz < < 1.1 THz in the 
small system has to be "reshuffled" to smaller frequencies leading to a broadening 
and an increase of the quasielastic line around u = 0. Since the quasielastic line is 
outside the frequency resolution of our Fourier transformation the consequences in 
the change of the quasielastic line can be observed better in the Fourier transform 
of Ssiq,!^), i.e. the incoherent intermediate scattering function Fs{q,t), which is 
shown in Fig. |^ for the system sizes A^ = 8064 and A^ = 1008 at q = 1.7 A~^. 
We recognize from this figure that Fs{q,t) shows a two step relaxation behavior 
similar to the case of silica and fragile glassformers (Ngai, Riande, and Ingram 
1998). As expected from our results for Ss{q,i^) the scattering functions Fs{q,t) 
have no system size dependence for sodium. In Fs{q,t) for silicon and oxygen the 
height of the plateau increases and the a relaxation process shifts to longer times 
with decreasing system size. Furthermore, the scattering functions for the small 
system show a pronounced oscillation for t > 0.2 ps which is due to the fact that 
the boson peak excitations present in the small system cause a peak in Ssiqji^) 
whereas in the large system only a shoulder is present. Finally we mention that 
the finite size effects in the dynamics of SDS are found in the whole q range and, 
moreover, do not affect the static properties of SDS. 
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Figure 1: Comparison of the static structure factor from our simulation 

(solid line) with the experimental data of Misawa et al. (1980) (dashed line). 
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Figure 2: Distributions of the coordination number for SDS at T = 300 K (open 
circles) and T — 2100 K (filled circles), a) Po-Si(-2), b) PNa-o(-2) and PNa-Na(-2)- 
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Figure 3: a) Self part of the dynamic structure factor for silicon, sodium and oxygen 
for the system sizes N = 1008 and N = 8064 at the temperature T = 2100 K 
and the wave-vector q — 1.7 A~^. The vertical hne at u — 1.1 THz marks the 
frequency below which there is a loss of intensity in the small system, b) Incoherent 
intermediate scattering function Fs{q,t) under the same conditions as in a). 
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